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ABSTRACT 

We develop a framework for studying the statistical properties of current sheets in numerical sim- 
ulations of 3D magnetohydrodynamic (MHD) turbulence. We describe an algorithm that identifies 
current sheets in a simulation snapshot and then determines their geometrical properties (including 
length, width, and thickness) and intensities (peak current density and total energy dissipation rate). 
We then apply this procedure to simulations of reduced MHD turbulence and perform a statistical 
analysis on the obtained population of current sheets. We evaluate the role of rcconnection by sepa- 
rately studying the populations of current sheets which contain magnetic X-points and those which 
do not. We find that the statistical properties of the two populations are different in general. We 
compare the scaling of these properties to phenomenological predictions obtained for the inertial range 
of MHD turbulence. Finally, we test whether the reconnecting current sheets are consistent with the 
Sweet-Parker model. 

1. INTRODUCTION 

Magnetic reconnection and turbulence are two of the 
most important and well studied processes in modern 
plasma physics and plasma astrophysics. Both of them 
have been recognized as critical to understanding many 
laboratory, space, and astrophysical phenomena and 
have been the subjects of extensive research in the past 
few decades. However, in many real natural systems, 
these processes do not occur in isolation, independent of 
each other. Real systems are often complex; they usu- 
ally involve many interacting components and many pro- 
cesses acting simultaneously. It is therefore clear that, 
in order to understand how reconnection and turbulence 
operate in such systems, one needs to understand how 
they interact and affect each other. 

In principle, one can define two sets of questions: (1) 
how a relatively small-scale externally superimposed or 
internally-generated turbulence affects reconnection of a 
large-scale magnetic field; and (2) how reconnection at 
the bottom of the turbulent cascade affects/controls en- 
ergy dissipation in magnetohydrodynamic (MHD) turbu- 
lence. Despite the obvious importance of both of these 
fundamental physics questions, so far relatively little re- 
search has been done on them, compared with the num- 
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ber of studies dealing with either turbulence or reconnec- 
tion separately. 

On the first topic, we note that, even though the po- 
tential importance of externally imposed magnetohydro- 
dynamic turbu lence in reconnection has been recognized 
early on (e .g., 'Matthaeus &: Lamkin| 1986[ Lazarian & 



|Vish niac|pr999 Kim & Diamond 2001), most of the tradi- 
tional studies of reconnection have assumed a completely 
smooth and laminar background/initial state. Only rela- 
tively recently, with the advent of more powerful numer- 
ical simulations, has this topic received a p roper atten- 
tion (Kowal et al. 2009 Loureiro et al. 2009). In addition 
to these studies, which focused on the ettect of an exter- 
nally driven turbulence on reconnection, one can also ask 
whether and how a reconnection process is affected by 
internally-generated MHD turbulence, produced as a by 



produ ct of the reconnection process itself (e.g., Strauss 



19881. This question has also received a lot ot atten- 



tion m recent years, largely in connection with the tran- 
sition of a new dynamic, non-stationary regime of fast 
resistive-MHD reconnection mediated by the seco ndary 
tearing/plasmoid instability ( Loureiro et al.| 2007|. This 
so-called plasmoid reconnection regime has been studied 



both analytically ( L oureiro et al. 2007 Uzdensky et al. 
2010| iBaalrud e t al. 2012: Loure iro et al.||2013j) and nu- 
merically (|Lapenta 2008; Samtaney et al.||2009| |Bhat- 



tacharjce et al. 200 9t ICassak et al . 2009; Huang & BhaF 
tacharjeen2010( |Loureiro et al.||2012^ and can be seen as 
to 



a special form ot turbulent reconnection. 
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As for the second topic — turbulent energy dissipation 
via magnetic reconnection at the bottom of the turbu- 
lent cascade — the number of papers published on it so 
far has been much less. The most notable was a recent 
series of statistic al studies by Servidio et al. (Servidio 
et al.|[2009, 2010l |2011bp] ). There has also been a more 



general investigation of energy dissipation in intermit- 



tent current and vorticity structures performed by Urit- 
sky et al. (2010). We discuss these papers in more detail 



and compare them with our analysis in the last section 
of this paper. 

In addition to the obvious fundamental-physics moti- 
vation for investigating reconnection in turbulence, there 
are also practical applications to various astrophysical 
environments. At issue here is the larger question of the 
intermittency of energy dissipation and hence of plasma 
heating in MHD turbulence. This question becomes es- 
pecially important in situations where there is strong 
prompt radiative cooling that may evacuate the dissi- 
pated energy from localized dissipation sites (hot spots) 
faster than it can be redistributed uniformly across the 
plasma by turbulent diffusion or thermal conduction. 
This can greatly alter the thermodynamics of such sys- 
tems; in particular, lead to a highly inhomogeneous tem- 
perature structure. An important astrophysical exam- 
ple of a system where this may happen is coronae of 
accretion disks around black holes accreting at a near- 
Eddington limit, e.g., in quasars, where ~ 100 keV elec- 
tron s are subject to very pow erful inverse-Compton cool- 
ing ( Goodman k. Uzdensky|2008j . However, these issues 
may also be important in various other high-energy as- 
trophysical systems, such as accretion disks and flows, 
hot X-ray gas in some galaxy clusters, etc. Furthermore, 
even in low-energy-density astrophysical systems, where 
prompt radiative cooling is not important, the issue of 
intermittency of turbulent energy dissipation may still 
be of interest for another reason. Namely, it determines 
the spatial distribution and coherence lengths of electric 
fields, and thus may affect the effectiveness of nonther- 
mal particle acceleration in coUisionless or weakly col- 
lisional plasmas in, e.g., radiatively-inefhcient accretion 
flows, galaxy clusters, and the solar wind. 

Another application is the dissipation and intermit- 
tency of the solar wind. Cur rent sheets are ass ociated 
with magnetic discontinuities ( Greco et al.||2010 ), which 
are observed in-situ by spacecraft in the solar wind and 
have been studied in a number of works recently. The 
origin of the discontinuities is not well understood and 
attracts considerable interest. Two possible explanations 
are discussed in the literature. One possibility is that the 
discontinuities are generated in_the solar corona an d then 



advected by the solar wind ((Borovsky 2008 ; Miao et al 
2011 1 . The other possibility is that they are dynamically 
generated in the solar wind, fo r example, due to inherent 



magnetic plasma turbulenc e (Li et al. 2011 Borovsky 
Zhdankin et al. 2012). I'he investigation of cur- 
rent sneets formed in turbulence, and their relationship 
to observable magnetic discontinuities, may therefore be 
useful for distinguishing between these two possibilities. 

The present paper is devoted to a detailed study of the 
intermittency of magnetic energy dissipation in magneto- 
hydrodynamic (MHD) turbulence and, more specifically, 
to assessing the role of small-scale magnetic reconnec- 
tion events at the bottom of the turbulent cascade in the 



overall turbulent heating process. For simplicity, in this 
study we focus only on small-scale magnetic structures 
and associated ohmic heating, leaving viscous dissipation 
of kinetic energy to a future study. Since ohmic heating 
per unit volume is equal rjj^ and we assume constant re- 
sistivity ?7, the main sites of magnetic energy dissipation 
are the regions of concentration of the current density j. 
In three-dimensional MHD turbulence such regions are 
two-dimensional current sheets, which may or may not 
be associated with magnetic reconnection (i.e., with a 
change of magnetic field topology, marked by the pres- 
ence of magnetic X-points or X-lines) . 

Correspondingly, the main focus of the present paper is 
on a statistical study of intense current sheets and their 
associated heating rates in two- and three-dimensional 
(2D and 3D) MHD turbulence. In this study we are 
interested in looking at current sheets as a population 
and want to address statistical questions such as: How 
important are intense current sheets in the overall en- 
ergy dissipation? What is the distribution of the current 
sheets in intensity (i.e., the peak values of the volumet- 
ric current density, jmax) and in integrated energy dissi- 
pation/heating rates? What is the distribution of their 
geometric properties such as lengths, widths, and thick- 
nesses? How are all these quantities correlated with each 
other? Are there substantial statistical differences be- 
tween current sheets that are associated with reconnec- 
tion events (containing an X-point) and those that are 
not? What is the role of discrete reconnection events in 
the overall magnetic energy dissipation and in its inter- 
mittency properties? 

To address these questions, we develop a set of numer- 
ical procedures and tools for identifying current sheets in 
simulation data and quantiatively characterizing them in 
terms of their peak current density, dimensions, overall 
magnetic dissipation rate, etc. Our algorithm allows us 
to distinguish between reconnecting current sheets (those 
with X-points) and non-reconnecting ones (without X- 
points). We then apply these tools to post-process our 
existing high-resolution (up to 1024^ grid points) numer- 
ical simulations of driven MHD turbulence. 

As a result of this study, we are able to show that a 
large number of the current sheets do not contain recon- 
nection sites, and likewise, many reconnection sites do 
not reside in intense current sheets. However, the cur- 
rent sheets that do contain reconnection sites tend to be 
larger and more intense, and their properties show more 
robust scaling relationships than current sheets without 
reconnection sites. We find that the scalings differ from 
what is expected from phenomcnological modeling of the 
inertial range of MHD turbulence, apparently refiecting 
the fact that current sheet formation also depends on 
the processes in the dissipation range of the energy spec- 
trum. A comparison for reconnecting current sheets of 
measured thickness with the thickness expected from the 
Sweet-Parker model shows reasonable agreement. 

This paper is organized as follows. In Section [2] we 
discuss phenomcnological estimates for the dissipation 
scale in MHD turbulence. In Se ction [3j we describe our 
numerical MHD simulations (§ |3.1[ ) and the algorithm 
we developed to identify and characterize current sheets 
(§ 3.2). In Section 111 we describe the results of our sta- 
tistical analysis of tne current sheets and compare to the 
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Sweet-Parker theory. In Section [5] we compare our re- 
sults to similar studies, discuss the potential implications 
for various space and astrophysical systems, and outline 
future extensions of this work. Finally, in Section [6] we 
summarize our conclusions. 

2. FIDUCIAL DISSIPATION SCALE OF MHD TURBULENCE 
The incompressible MHD equations take the form 



!=F . V) z± = -VP + i^V^z* + f=^, 



V-z±=0 (1) 

where = v ± b are the Elsasser variables, v is the 
fluctuating plasma velocity, b is the fluctuating magnetic 
field (in units of the Alfven velocity, Va = Bo/-^47r/9o, 
based upon the uniform background magnetic field Bq 
and plasma density po)j P = {p/po + &^/2), P is the 
plasma pressure, p is the fluid viscosity (which, for sim- 
plicity, we have taken to be equal to the magnetic dif- 
fusivity), and represent forces that drive the tur- 
bulence at large scales. It can be shown that in the 
limit of small amplitude fluctuations, and in the ab- 
sence of forcing and dissipation, the system describes 
non-interacting linear Alfven waves with the dispersion 
relation w^(k) = ifcyVyi. The incompressibility condi- 
tion requires that these waves be transverse. Typically 
they are decomposed into shear Alfven waves (with po- 
larizations perpendicular to both Bq and to the wave- 
vector k) and pseudo-Alfven waves (with polarizations 
in the plane of Bq and k and perpendicular to k). 

In order to set the stage for the subsequent sta- 
tistical study of intense intermittent dissipative struc- 
tures in MHD turbulence, in this section we discuss 
the phenomenological estimates for the dissipation scale 
and, more generally, for the typical turbulent dissipative 
structures. These estimates are based on the analysis 
of the scaling relationships presented in ( Boldyrev||2005 



20061 that characterize the inertial range of the turbu- 



lence and on evaluating the conditions when dissipative 
(e.g., resistive) terms become important. Importantly, 
these estimates are done without taking into account in- 
termittency of MHD turbulence and thus can serve as 
the simplest baseline against which intermittent dissipa- 
tive structures, which constitute the main focus of this 
paper, can be compared. 

According to the ph enomenological picture of M HD 
turbulence discussed b y Goldreich & Sridhar ( 1995 ) and 
Boldyrev ( 2005 2006 ) , the typical structures within the 
inertial range are highly 3D-anisotropic. The degree of 
anisotropy of a turbulent eddy [which is related to the 
scales over which the typical (rnis) fluctuations of veloc- 
ity and magnetic field are correlated in different direc- 
tions] increases as o ne goes down to smalle r and smaller 
scales. As argued by Boldyrev (2005 20061, it is charac- 
terized by the following relationships between the typical 
eddy scales in different directions: 



(2) 
(3) 

Here, I is the scale of a given structure in the field- 
parallel direction, ^ is the typical scale in the field- 
perpendicular direction along the fluctuating velocity 



and magnetic fields and b^, and finally, A is the typ- 
ical scale in the field-perpendicular direction but across 
Va and Thus, we see that indeed, the turbulence is 
highly anisotropic, I ^ ^ X. 
The typical velocity an d magnetic field fluctua tions in 



the inertial range scale as Boldyrev ( 2005 2006 ) 



vx (X bx (X 



(4) 



Furthermore, from this we can obtain scalings for sev- 
eral other key quantities, e.g., the characteristic eddy 
turn-over time (this can be obtained from Goldreich- 
Shridhar critical balance argument): 

Tx (X Ix/Va cx Cx/ix cx A^/^ , (5) 
and the characteristic current density at a given scale: 



jx (X — (X X . 



(6) 



It should be noted that phenomenological models typ- 
ically deal with an idealized picture which addresses the 
scaling of the fluctuating flelds in the limit of very large 
Reynolds and magnetic Reynolds numbers. In reality 
these numbers are limited, and they are relatively small 
in numerical simulation s (on the order of 10^). As was 
shown in (Wang et al.„201l| [Boldyrev et al.|2012b|a ) 



the case ot limited inertial mterval the scalmgs of mag- 
netic and velocity fluctuations appear to be slightly dif- 
ferent due to residual energy generated at large scales. 
The same effects ar e also observed in the solar wind tur- 
bulence (e.g., Boldyrev et al. 2011). In our present dis- 
cussion we do not take into account such effects, for two 
reasons. First, the phenomenological theory is not devel- 
oped well enough to address the scales near or inside the 
dissipation interval with the same certainty with which it 
addresses the inertial interval, and second, the statistical 
relations measured in the present work are not neces- 
sarily related to the second-order correlation functions 
predicted by most phenomenological models. 

Now let us turn to the discussion of the dissipation 
scale, which we will denote by A,;. Similarly, all the quan- 
tities evaluated at the dissipation scale will be denoted 
by the subscript 77 (assuming for definiteness/simplicity 
that viscosity ly is equal to the magnetic diffusivity 77, i.e., 
that the magnetic Prandtl number is 1). The scale A^ is 
defined as the scale at which the resistive diffusion time 
'tj — ^Iq/'n^ is comparable to the corresponding 



across it, t. 



eddy turn-over time, which gives 

A, cx 7^2/3 , (7) 

and, correspondingly, 

a r,^/'^ , cx 
The other turbulent quantities at this scale are estimated 

^ Note that the relationship I oc ^^/-^ would correspond to 
the Goldreich-Shridhar (1995) theory, although the scaling of the 
fluctuatin g field 6^ Ei^nd vx is different in the phenomenology of 
|Boldyrev||2006|l , which we consider here and which agrees better 
witn tne numerical data. 
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Reduced MHD (RMHD) model: 
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'1/2 



1/6 



(9) 

(10) 
(11) 



Finally, the characteristic energy dissipation rate in 
a single typical dissipative structure can be estimated 
(in 3D) as 

£^ ~ 77j2 X^^^l^ (X 77^/2 . (12) 

For reference, since most current state-of-the-art 3D 
numerical simulations of MHD turbulence have ry ~ 10"'^ 
(in units normalized to the energy-containing length- 
scale and the rms velocity/magnetic field at that scale), 
the typical values of the above fiducial dissipation-scale 
parameters are: A,, ~ 0.01, ~ 0.03, ~ 0.1, 
^ 6r, ^ 0.3, T,j ~ 0.1, jrj ^ 30, and ^ 3 x 10~^. 

Next, it is interesting to note that the above estimates 
are consistent with the notion that such typical dissi- 
pative s tructures can be vi ewed as elementary Sweet- 
Parker (Sweet 1958 Parker|,1957j reconnection current 



sheets, m the sense that their lifetime is compara- 
ble to the Sweet-Parker (SP) reconnection time for the 

amount of flux equal to A,, 6,,, as one can easily see. The 
characteristic Lundquist number of these current sheets 

is 5*^ oc brj^ri/v cx r]~^^^ cx where 5*0 oc 7y~^ is the 

global Lundquist number. Since in most modern sim- 
ulations of 3D MHD turbulence this global Lundquist 
number does not exceed about 10,000, the Lundquist 
number of the typical SP current sheets at the dissi- 
pative scale is typically of order Sr/ ~ 20 or smaller, 
and correspondingly, their aspect ratio is not very large, 

5*^^ ~ Sg^^ ~ 4. Thus, it is not surprising that such 
current sheets are usually not recognized as thin SP cur- 
rent sheets in the simulation data. Their character as 
SP current sheets may have been more clearly visible if 
we had access to simulations with much larger (not 
possible in the near future). 

The estimates obtained in this section represent 
the phenomenological predictions for typical dissipation 
structures in 3D MHD turbulence. As we shall see in 
this paper, while they may be responsible for a signifi- 
cant fraction of the actual dissipation, they do not ac- 
count for all of it; a sizable fraction of turbulent energy 
is dissipated in a few much more intense structures, and 
the main goal of this paper is to assess this contribution 
quantitatively. 

3. NUMERICAL PROCEDURES 
3.L Numerical MHD Simulations 



For strong MHD turbulence, Goldreich & Sridhar (Gol- 
dreich fc Sridhar] |T995) argued that the pseudo-Altven 
modes are dynamically irrelevant for the turbulent cas- 
cade (since strong MHD turbulence is dominated by fluc- 
tuations with k± ^ fcj|, the polarization of the pseudo- 
Alfven fluctuations is almost parallel to the guide field 
and they are therefore coupled only to field-parallel gra- 
dients, which are small since /cy <C fcj^). If one filters 
out the pseudo-Alfven modes by setting zjj' = 0, it can 
be shown that the resulting system is equivalent to the 



We note that in RMHD the fluctuating fields have only 
two vector components, but that each depends on all 
three spatial coordinates. Moreover, because the are 
assumed incompressible ( V • z^ = 0) , each field has only 
one degree of freedom; this is more commonly expressed 
in terms of stream functions in the more standard form of 



the RMHD e quations ( ,Kadomtsev fc Kantorovich 1974 
Strauss|1976 ). This form is obtained if we introduce the 



axial component of the vector potential = ~-^z E^nd the 
stream function (j) = xl Bq^ where x is the scalar poten- 
tial. The magnetic field and velocity are then recovered 
from and cf) via 



B = X 'S/lp + BqBz 

The reduced MHD equations are then written as: 



(14) 



P 



dip 

In 

duj 
dt 



pu ■ Voj 



Mo 



= 



1 



piyV^Lj = — B ■ Vj 



(15) 



where uj — V'^^ and j ~ V^V'- 

We solve the RMHD equations ( 13 ) in a periodic, rect- 
angular domain with dimensions L^x Ly , where the sub- 
scripts denote the directions perpendicular and parallel 
to the background guide field Bq, respectively. We set 
L± — 27r, L\\/L± — 6 and Bq = hhrms^-z, where forms 
is the root-mean-square average of the fiuctuating mag- 
netic field component. The turbulence is driven at the 
largest scales by colliding Alfven modes We drive both 
Elsasser populations by applying statistically indepen- 
dent random forces f and f ~ in Fourier space at wave- 
numbers 27r/L_L < fc_L < 2(27r/L_L), k\\ = 2t:/L\\. The 
forces have no component along z and are solenoidal in 
the xy-plane. All of the Fourier coefficients of out- 
side the above range of wave-numbers are zero and inside 
that range are Gaussian random numbers with ampli- 
tudes chosen so that Vrms ~ 1- The individual random 
values are refreshed independently on average approxi- 
mately 10 times per turnover of the large-scale eddies. 
The variances a\. — (|f*P) control the average rates 
of energy injection into the z+ and fields. In this 
work we consider the "balanced" case, that is we choose 
a'^ = u~ . 

A fully dealiased 3D pseudo-spectral algorithm is used 
to perform the spatial discretization on a grid with a 
resolution of A'^ x A^n mesh points (we typically take 
= N\\ = 512 or 1024, see Table [l]). We note that the 
domain is elongated in the direction of the guide field in 
order to accommodate the elongated wave-packets and 
to enable us to drive the turbulence in the strong regime 
while maintaining an inertial range that is as extended 

^ Turbulence can also be driven by driving v or b fluct uations 
at large sca les: this does not affect the inertial interval, see | |Mason| 



at large sca li 
|et al. [20081 . 



(13) 
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as possible (see Perez fc Boldyrev]|2010 ) . This is a phys- 
ical requirement that should be satistied no matter what 
model system, full MHD or reduced MHD, is used for 
simulations. We also note that magnetic Reynolds num- 
ber is defined by i?^ — b^,^s{L^/2'K) /v in our studyj^ 

3.2. Numerical Procedures for Identifying and 
Characterizing Current Sheets 

The question that wc consider is: how can we unam- 
biguously identify current sheets in a numerical simula- 
tion, where the thinness of developing current sheets is 
always limited by the grid resolution? To address this 
issue, we develop a concrete methodology and a specific 
algorithm for identifying structures in a large 3D simula- 
tion and for characterizing them quantitatively in terms 
of their geometrical properties and intensities. 

Our algorithm is designed to be impartial to reconnec- 
tion is the sense that it does not automatically associate 
the dissipation sites with magnetic X-points. Hence, the 
current sheets may or may not coincide with magnetic 
X-points. We consider the relationship to magnetic re- 
connection only after the current sheets have been fully 
characterized b y their other propert ies. This differs from 
the approach of Servidio et al. ( 2010 1, where reconnection 
sites are detected first and the corresponding dissipation 
regions are studied afterwards. 

We implement the algorithm in IDL (Interactive Data 
Language). The algorithm can be applied either to 2D 
simulations or to 3D simulations with a strong magnetic 
guide field Bq in the z-direction; here, we focus on the 
3D case since it is more general. We assume that current 
sheets are predominantly aligned with the background 
guide magnetic field in order to simplify some parts of 
the algorithm (this assumption is supported by the re- 
sults). The required input quantities are the magnetic 
flux function ijj, the current density j = jz = V^V; 
and the perpendicular {xy) magnetic field b — z x Vif). 
An illustration of the current density profile in an xy- 
plane cross section of our data is shown in Fig. [T] where 
the presence of numerous current sheets penetrating the 
cross section is evident. 

Before applying the algorithm, it is useful to inter- 
polate the data to a higher-resolution grid by using a 
Fourier zero-padding and interpolation technique ( |Ser-| 
vidio et al. 2010). This procedure consists of applying a 
fast i<burierTransform (FFT) on the data, zero-padding 
the resulting set of coefficients to produce a larger array, 
and applying an inverse FFT to obtain the original points 
along with interpolated points in between. This interpo- 
lation scheme cannot be used on an entire 3D dataset 
because the computational costs are too great. However, 
it can be applied on any 2D cross-section of the data, 
from which we obtain a AN^ x ANy array in place of 
the original x Ny array. This is useful because it im- 
proves the accuracy of several of our measurements, such 

^ In the case of reduced MHD though, when the components 

are exphcitly removed, the resulting system l |13[ l is invariant with 
respect to simultaneous rescaling of the background field Bq and 
the field-parallel spatial dimension of the system, if one neglects the 
dissipation terms. Therefore, one can rescale the field-parallel box 
size to L|| = L^, that is, conduct the simulations in a cubic box, 
provided the backgrond field Bq is rescaled a ccor dingly. We should 
note however that the dissipation terms in are not invariant 
and they should be changed accordingly under such rescaling. 
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Figure 1. The current density in an a;j;-plane cross section of 
data. Red indicates negative current and blue indicates positive 
current. Current sheets protruding through the cross-section are 
clearly present. 

as width and thickness, which are made in xy planes. 

We now describe our current-sheet identification algo- 
rithm. The first step is to locate the current sheets in a 
given simulation time snapshot. Since current sheets are 
characterized by extrema in the current density profile, 
the problem reduces to finding local maxima in the (mag- 
nitude of) current density. To achieve this, the algorithm 
scans through all points above a pre-specified threshold 
current density jthr (which is several times larger than 
the globally averaged magnitude of current density) and 
selects those points that are local maxima within a sur- 
rounding cubic subarray of the data. Each subarray is 
centered at the candidate point and has size (2n -|- 1)^, 
where n is a parameter. With larger n, the algorithm 
finds only the most dominant peaks. We choose n = 4 
for our analysis, which allows well-resolved peaks to be 
detected. Our results do not change significantly if n is 
adjusted by a factor of two or so. Every maximum is 
then identified with a current sheet that we label by in- 
dex I, and the corresponding current density is referred 
to as the peak current density, jmax,i- 

The second step is to identify the points belonging to 
each given current sheet. These are defined to be the 
points that collectively connect to the point of peak cur- 
rent density, with the condition that each point has a 
magnitude of current density greater than a minimal 
value, jmin,i- For definiteness, we choose the current 
sheet boundary to be half of the peak current density, 
so jinin,i = Jmax,i/2. Thc algorithm determines the cur- 
rent sheet points in the following way. First, it considers 
the points adjacent to the peak (four points for 2D anal- 
ysis, six points for 3D analysis), and from these points 
adds the ones with current densities greater than jmin.i 
to a list. Then for each point now on the list, the al- 
gorithm adds to the list any unchecked adjacent points 
with current density above jmin,i- This procedure is re- 
peated for subsequent points added to the list until no 
new points are found. Provided that jthr is chosen to be 
high enough, the current sheets thus constructed are rel- 
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Figure 2. Highlighted in green are the xy-plane cross sections 
of current sheets found by applying the algorithm to 3D current 
density data. The cross section of data is the same one as used in 
Fig.Il] Note that the dominant structure is that of a current sheet. 

atively isolated and sparse and do not form one globally 
percolating cluster. 

There is, however, an ambiguity regarding how to treat 
current sheets that contain multiple peaks. If a partic- 
ular current sheet has jmin.i > Jthn then it is possible 
for a nearby peak to be associated with a second current 
sheet that contains points shared with the first current 
sheet. Whether or not such current sheets should be 
regarded as independent is unclear - this ambiguity is 
similar to the to pography prob lem of objectively defin- 
ing a mountain (Gerrard 1990). We choose to redefine 
current sheets with overlapping boundaries (by the crite- 
ria above) to be a single current sheet, with jmax.i taken 
from the most intense peak and jmin,i taken to be half 
of the smallest discernable peak within the composite 
current sheet. Therefore, a composite current sheet will 
have jmin,i < Jmax,i/2. In general, the statistical prop- 
erties and scalings of current sheets should be largely 
independent of how exactly the composite current sheets 
are defined. 

An illustration of the resulting cross sections for cur- 
rent sheets in our data is shown in Fig. [2j Additionally, 
the current sheet cross sections in the xz plane (parallel 
to the guide field) for 1 /3 of the simulation box height are 
shown in Fig. [Sj It is clear that the dominant structure 
is indeed that of a current sheet — a quasi-lD structure 
in a 2D slice or a quasi-2D ribbon in full 3D. Upon closer 
inspection, the actual 3D shape of a current sheet can 
often be complex. Common structural features include 
curvature, irregular boundaries, and strong asymmetry 
around the peak. This departure from the traditional, 
idealized Sweet-Parker picture of a straight and smooth 
current sheet can be attributed to turbulence. To some 
degree, the structure of an individual current sheet de- 
pends on the criteria in our definition; a current sheet 
that looks irregular when using the algorithm with one 
set of criteria may appear more regular when applying 
another set of criteria. However, the statistical conclu- 
sions of our study are not very sensitive to these criteria. 
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Figure 3. On the left, an x2-plane cross section of the current 
density is shown, for 1/3 of the simulation box height with actual 
aspect ratio. On the right, areas corresponding to current sheets, 
as found by applying the algorithm to the 3D current density data, 
are highlighted in green. Note that the current sheets are elongated 
in the direction parallel to the background guide field. 



With each individual current sheet now associated with 
a set of points in space, we can characterize their main 
physical properties quantitatively. In the following, we 
describe how we compute the ohmic energy dissipation, 
width, thickness, length, and upstream magnetic field 
vectors for a given current sheet. Note that we order the 
three current sheet dimensions by A < ^ < /, where A is 
thickness, ^ is width, and I is length. 

We start with f , the total magnetic energy dissipated 
by the current sheet per unit time. Note that we con- 
sider only ohmic dissipation while ignoring viscous dis- 
sipation, which, in principle, can be of comparable mag- 
nitude in our simulation since we have Pm = v jr] = 1. 
Our diagnostic procedures can be generalized to include 
viscous dissipation rather straightforwardly and this is 
something we are planning to do in the near future. Since 
integration over the current sheet volume is numerically 
equivalent to a summation over all N points, the energy 
dissipated becomes 



VJ' 



N 
k=l 



(16) 



hy, and hz are the dimensions of 



where jk is the current density at the fcth point of the 
current sheet, and h 
each cell. 

The measurements of thickness and width for each cur- 
rent sheet are made on the xy plane that contains the 
point of peak current density of that current sheet. Ac- 
cordingly, it is useful to define the current sheet cross 
section to be the set of points in this plane that belong 
to the current sheet. 

A relatively direct method for measuring current sheet 
thickness is as follows. We first determine the direction of 
most rapid descent from the peak current density. This 
is accomplished by numerically computing the Hessian 
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matrix for the current density j{x,y) at the peak, 



dyxj dyyj 



(17) 



and calculating the eigenvectors of H. We then find the 
distance from the peak in this direction at which the 
current density drops to that of the boundary, jmin,i- We 
repeat this procedure in the opposite direction, and add 
the two distances to get the total sheet thickness A. We 
also employ an alternative way to estimate the thickness 
(similar to the one used by Uritsky et al. ( 2010 )) in which 
we divide the area of the current sheet cross section by 
its width (for which the measurement procedure is dis- 
cussed in the next paragraph). A scatter plot compar- 
ing the thicknesses obtained by using the two different 
methods is shown in Fig. [4] for simulations with mag- 
netic Reynolds number — 1800 and Rm — 3200 with 
resolution 1024^. In general, we find a good agreement 
between these two independent methods. The increased 
amount of scatter for the Rm = 3200 case, which has 
thinner current sheets than Rm = 1800, suggests that 
those current sheets may not be sufficiently well resolved. 
Any measurement of thickness is intrinsically limited by 
the numerical resolution of the simulation. Indeed, the 
thickness obtained from the eigenvector method is quan- 
tized due to stepping across an integer number of cells of 
width h (or h/A when interpolation is used). The other 
thickness estimate is less affected by resolution, but loses 
information on the local thickness near the peak, the lo- 
cation of most interest. 

Next, we define the width ^ to be the largest distance 
between any pair of points in the xy cross section of the 
current sheet — this is very accurate unless the cross 
section is highly curved. It is found by iterating over all 
pairs of points in the cross section to obtain the maxi- 
mum distance. It is worth noting that the width of the 
current sheet can be found alternatively by applying the 
eigenvector measurement procedure used for thickness, 
but using the other eigenvector of H instead. However, 
this method is less robust since the boundary can be 
reached prematurely due to current sheet curvature. 

One may consider defining the length I of the current 
sheet to be the maximum distance between any pair of 
points in the entire 3D current sheet ribbon. However, 
this is not useful in practice because it requires iterating 
over all pairs of points in the current sheet, which is 
computationally expensive for large 3D current sheets. 
Instead, we use the assumption that current sheets are 
aligned in the z-direction, and define the length to be the 
distance between the endpoints, which are defined to be 
the current sheet points with the maximum z coordinate 
and the minimum z coordinate. This estimate is accurate 
even for current sheets that are misaligned with respect 
to the z-direction, because we take into account the x 
and y coordinates of the endpoints. 

We also determine the magnetic field vector at key 
points of the current sheet. In particular, we measure the 
asymptotic (upstream) values of the xy magnetic field on 
both sides of the current sheet. We define Si y and i?2.j| 
to be the magnetic field components parallel to the cur- 
rent sheet at the two edges, which will be used later to 
compare to reconnection models. We first measure the 
magnetic field vectors Bi and B2 at the two edge points 



which we define to be the points used in the thickness 
measurement (see above). Then, and B2^\\ are found 
by projecting Bi and B2 onto the direction perpendic- 
ular to the eigenvector of H in the xy plane used to find 
the thickness. If v is the unit vector in the xy plane that 
is orthogonal to the eigenvector of H, then 

=Bi V 

=B2 V (18) 

Finally, we consider the degree of association of the 
current sheets with magnetic reconnection events. For 
the purposes of this paper, we use magnetic X-points 
as a proxy for reconnection sites, although 3D reconnec- 
tion is not alwa ys associated with n ull points (Priest & 
Demou lin IQQSj' jParnell et al.|[20Tol ). Thus, we first de- 
tect X-points in the simulation, and then classify current 
sheets by whether or not they contain an X-point. 

We note that X-points are equivalent to saddle points 
in the magnetic flux function. Therefore, the problem 
in 2D reduces to detecting saddle points and determin- 
ing whether they lie inside of current sheets. The 3D 
case is subtler because there is no analogue of a saddle 
point for a function in 3D space. However, we can again 
take advantage of the fact that reduced MHD is used, 
which allows us to consider saddle points in the mag- 
netic flux function for every xy-plane that constitutes 
the volume. The set of these X-points then approximates 
X-lines where 3D reconnection takes place. 

One way to find saddle points is by using the first and 
second derivative test. In this case, saddle points corre- 
spond to points where the two eigenvalues of the Hessian 
matrix for ijj{x,y), H^p, have opposite signs. We found, 
however, that this method did not work well for discrete 
numerical data. As a simple example, take a region of 
space where the function has integer values: 

2 2 
2 1 
2 

The central point should be considered a saddle point 
because the function alternatively decreases in two direc- 
tions and increases in two orthogonal directions. How- 
ever, the slope does not change sign in the x-direction 
nor the y-direction, and so this saddle point would be 
missed by the first derivative test. 

We use an alternative method that always finds these 
saddle po ints, similar to one proposed by Kuijper (Kui- 
jper|2004 ). The algorithm first considers the eight points 
surrounding the candidate point. If the function at 
these points, in clockwise order, rises above the candi- 
date point's value twice and falls below the candidate 
point's value twice, then it follows that the candidate 
point is a saddle point. If the function rises and falls any 
other number of times, then it is not a saddle point. 

A disadvantage of this procedure is that some false sad- 
dle points may be found if there is small-scale structure 
near the point. To avoid this, we repeat the procedure for 
two larger clockwise loops, having circumferences of six- 
teen and twenty- four points, respectively. This is shown 
schematically in Fig. [51 If the pattern of rising twice and 
falling twice still holds for each of these loops, then we 
consider the candidate point to be a saddle point. After 
determining the location of all of the saddle points, it is 
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Figure 4. A scatter plot comparing the thickness measured directly by the algorithm (on the j/-axis) and the thickness estimated by 
dividing the current sheet cross sectional area by width (on the x-axis). The left panel shows the simulation with Rm = 1800 while the right 
panel shows the simulation with i?^ = 3200. The significant amount of scatter between the two methods demonstrates the importance of 
choosing an appropriate method. Also, the increased amount of scatter in the Rm, = 3200 case may suggest that the current sheets are not 
sufhciently resolved. 

jjj ~ 30. At the same time, we always choose jthr to be 
significantly smaller than the global simulation box max- 
imum, max(j). For any given simulation, the lower limit 
for jthr is set by the value at which connected islands 
of high current density span the entire simulation box 
(similar to problems in percolation theory) . By choosing 
the values of jthr shown in Table [TJ the largest structures 
have a length comparable to the dox size in the z direc- 
tion. Also shown in Table [T] is the average number of 
current sheets detected per snapshot. The total number 
of current sheets in the sample are Ntot — 7175 for Case 

Figure 5. To determine whether a candidate point is a sad- 1^ Ntot = 6616 for Case 2, and Ntot = 9343 for Case 3. 

die point, consider the surrounding points along a loop and check 
whether those values twice rise above and twice fall below the value 

at the candidate point. To avoid detecting false saddle points, re- rr, v,i 

peat the procedure for two larger loops shown. laDle 

simulation parameters 

straightforward to classify the current sheets by whether 
or not they contain an X-point, and are therefore asso- 
ciated with magnetic reconnection. 

This concludes the discussion of the algorithm. We 
have identified the locations and volumes of current 
sheets, measured several of their properties, determined 
the nearby magnetic field vectors, and classified them by 
whether or not they contain an X-point. The next sec- 
tion discusses the results that we obtain from applying 
the algorithm to our simulation data. 

4. RESULTS 

We analyze three different simulations of MHD turbu- 
lence. In Case 1, we use resolution of 512'^ and magnetic 
Reynolds number — 1800. In Case 2, we use resolu- 
tion 1024^ and Rm — 1800. In Case 3, we use resolution 
1024^ and Rm — 3200. We analyze multiple time snap- 
shots for each simulation (10 for Cases 1 and 2, and 7 
for Case 3). We find that the current sheets are best 
resolved in case 2, so we use data from that simulation 
for the majority of our analysis. The threshold current 
density for detection, jthr, differs for each shown 
in Table [T] We always choose jthr much larger than the 
mean, javg ~ 10, and than the fiducial estimate ( |11[ ) 
for the current density at the turbulent dissipation scale. 



Case 


Resolution 


Rm 


jthr 


max(j) 


A''/snapshot 


1 


512^ 


1800 


100 


437 


718 


2 


1024^ 


1800 


130 


656 


662 


3 


1024^ 


3200 


160 


837 


1335 



We find that many current sheets (more than half of 
the total) do not contain X-points and, likewise, many 
X-points do not lie in strong current sheets. This lack 
of correlation suggests two things. Firstly, turbulence 
creates current sheets that are not located at sites of re- 
connection. Secondly, active reconnection does not occur 
at many X-points in the domain. 

We find, however, that the strongest current sheets 
tend to have X-points, and that X-point containing 
sheets exhibit different statistics in general. It is there- 
fore reasonable to conclude that the strongest current 
sheets, and hence the most intense energy dissipation 
events are driven by reconnection, while the smaller ones 
are formed by turbulent fluctuations. 

We now present the details of our quantitative statis- 
tical analysis. Note that all of the following probability 
distributions are normalized such that the integral over 
the given values is equal to unity. Also, distribution plots 
are shown on a log-log scale to highlight power law re- 
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Figure 6. The probability distribution for current sheet ohmic 
energy dissipation rate, £, measured using our algorithm and 
shown on a log-log plot. The distribution for all current sheets 
(in black) shows a power law tail with an index near —1.8. Also 
shown are distributions for two subpopulations of current sheets: 
those with X-points (in red) and those without X-points (in blue, 
with a power law index near —2.7). Current sheets without X- 
points dominate at low energies, while current sheets with X-points 
dominate at high energies. 

gions. The spatial scales are measured with respect to 
the size of the simulation box in the direction perpendic- 
ular to the guide field. 

It is illuminating to first look at the mean values of 
the current sheet properties in the three simulations, 
even though the exact values depend on what detection 
threshold (jthr) is chosen. These properties are shown 
in Table [2| It is immediately seen that the structures 
are indeed highly anisotropic, i.e. (/) (^) (A). We 
also see that Case 1 and Case 2, which differ in resolution 
but have equal Reynolds numbers, exhibit similar length, 
width, and energy dissipation rates, despite having dif- 
ferent threshold current densities. We see however that 
the mean thickness, (A), changes significantly between 
the two cases, suggesting that it may be insufficiently 
resolved in Case 1. We also note that for all cases, the 
current sheets in total occupy less than 1% of the system 
volume but account for roughly 25% of the overall ohmic 
dissipation, showing that they are indeed a significant 
contribution to the total energy budget of the system. 



Table 2 

Mean values of current sheet parameters 



Case 


(0 


(0 


(A> 


( Jmax ) 




1 


0.39 


0.049 


0.0035 


130 


3.6 X 10-* 


2 


0.37 


0.046 


0.0024 


182 


3.5 X 10-* 


3 


0.29 


0.033 


0.0019 


213 


1.5 X 10"* 



4.1. Statistical analysis 

The probability distribution for energy dissipation rate 
integrated over a current sheet, £ from Eg. |16[ is shown 
in Fig. |6] The distribution for all current sEeets (black 
curve) exhibits a power law tail (for £" > 5 x 10""^) with 
an index of —1.8 ± 0.1. It is interesting to note that 
this power law is close to, but slightly harder than, the 
critical power law of index -2, indicating that the overall 
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Figure 7. The probability distribution for peak current den- 
sity, jmax, in current sheets. Current sheets with X-points (in red) 
dominate over sheets without X-points (in blue) at large current 
densities, in agreement with the asymmetry in energy dissipation 
rate (see Fig.lml. 

energy dissipation rate is dominated by large-f events. 
Also shown in Fig. [6] are similar probability distributions 
for the sub-populations of current sheets with X-points 
and current sheets without X-points (in red and blue, re- 
spectively). The current sheets with X-points dominate 
at large energy dissipation rates while the sheets without 
X-points dominate at small energy dissipations. In fact, 
every single current sheet with £ > 2 x 10^^ contains 
at least one X-point. Furthermore, approximately 84% 
of the total energy dissipation of current sheets occurs 
in sheets with X-points, while the remaining 16% occurs 
in sheets without X-points. The sheets without X-points 
exhibit a separate, steeper power law tail with an index 
close to —2.7. 

The probability distribution for peak current density, 
irnax, is shown in Fig. [7] In this case, there is a possible 
power law with index around 4, but is difficult to mea- 
sure because the range of current densities sampled is 
relatively narrow, less than an order of magnitude (due 
in part to our relatively high choice of current density 
threshold — we only have a factor of 5 between the 
threshold and the global maximum current density) . We 
find that the current sheets with large peak current den- 
sities tend to contain X-points, in agreement with the 
asymmetry in the distribution of £ for the two subpop- 
ulations. This supports the picture that the strongest 
sheets tend to contain X-points. 

The energy dissipation rate and peak current density 
are two different measures of a current sheet's strength, 
so we expect them to be correlated. Scatter plots of £ 
versus jmax are shown in Fig. [8 for the two populations of 
current sheets separately. We find that there is indeed a 
significant correlation between the £ and jmax for current 
sheets with X-points, with a power law dependence £ oc 
Jmax°'^- However, we find that the correlation for the 
current sheets without X-points is much weaker. 

Now we discuss the statistical properties of the current 
sheet sizes. The probability distributions for length I, 
width ^, and thickness A are shown in Fig.[9j For the pop- 
ulation of all current sheets, the distributions for length 
and width have power law tails, both with indices near 
—2.5. The thickness distribution declines much more 
rapidly at large A; if fit to a power law, this distribu- 
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Figure 8. As shown in the scatter plot on the left, there is a significant correlation between the ohmic energy dissipation rate and peak 
current density for current sheets that contain X-points, described approximately by a power law with index 3.5. However, as shown on 
the right, the correlation between the two quantities is much weaker for sheets that do not contain X-points. 
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Figure 9. The probability distributions for all current sheet dimensions. For the population of all current sheets (shown in black), the 
distributions for length / and width g have power law tails with indices near —2.5, while the thickness A has a much steeper tail (possibly 
an exponential, with a fit shown in green). The current sheets with X-points tend to be longer and wider (in red) than the sheets without 
X-points (in blue), but thickness is similarly distributed for all current sheets. 



tion "would have an index near —9, but the decHne is so 
steep that is not clear that a po"wer la"w is appropriate. 
An exponential decline appears to better match the data, 
"with a fit proportional to exp(— A/(3.6 x 10^^)) shown. 
The current sheets "with X-points tend to be longer and 
"wider than current sheets "without X-points, "while the 
distribution of thicknesses is similar for both types of 
current sheets. 

Next "we consider the cross-correlations o f th e current 
sheet dimensions. The scatter plots in Fig. [To]sho"w the 
correlation bet"ween all pairs of dimensions, for the pop- 
ulation of all current sheets. The length and "width are 
strongly correlated and have a nearly linear relationship, 
obeying a po"wer la"w I (x ^o.9±o.i_ This is some"what 
steeper than the index 2/3 naively expected from MHD 
turbulence (from Eq. [s]). This disagreement, ho"wever, 
is not surprising since the naive expectation describes 
typical turbulent eddies, and there is no a priori reason 
"why it should apply to the intense, strongly intermittent 
current sheets under investigation here. 

For the population of all current sheets, we find that 
thickness is not strongly correlated "with length or "width. 
There appears to be a broad trend for thickness to in- 
crease "with length and "width, but it is dominated by 
scatter. Note also that thickness measurements cover a 



shorter range of values (10"'^ ^ ^ 5 x 10~^) than the 
other t"wo dimensions, making it difficult to discern any 
trend. 

The picture of the scaling of dimensions changes some- 
"what "when "we consider separately the sub-populations 
of current sheets "with X-points and current sheets "with- 
out X-points. We find no significant difference bet"ween 
the t"wo populations "with regard to the scaling of length 
with width. However, the scaling with thickness does 
change. We find that both length and width are more 
strongly correlated with thickness for sheets that contain 
X-points, as shown in the left panels of Figs. 11 and [121 
The power-law fits for both of these correlations yield 
similar indices, I oc A'*"°^^"° and ^ oc X'^-^^^-^, Both of 
these scalings are much steeper than the scalings naively 
expected from MHD turbulence, which are I cx and 
^ C!c A'^/'* (see Eq.jsjand Eq.[2]). This substantial disagree- 
ment between observations and the phenomenologically 
expected power laws could be explained by the fact that 
thicknesses lie within the dissipation range, where the 
MHD turbulence estimates do not apply. For current 
sheets without X-points (see the right panels of Figs. 11 



and 12), we find that the correlations of lengths and 
widths with thickness are not as strong. A power-law 
fit gives I cx A'^"°*^"" and ^ cx A^"^*^"°, which is signifi- 
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Figure 10. For the population of all current sheets, the length I and width ^ are strongly correlated with a power law fit of index 
0.9 it 0.1 (in black). This is somewhat steeper than the naively expected index of 2/3 (in red). The current sheet thickness A is not strongly 
correlated with length or width. 
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Figure 11. The current sheet length I versus thickness A shown separately for sheets with X-points (left) and sheets without X-points 
(right). The sheets with X-points appear to have a stronger correlation between length and thickness, with a power law index near 4.0it 1.0. 
The length and thickness are not strongly correlated in sheets without X-points, with a possible power law index near 2.0 ± 1.0. 
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Figure 12. The current sheet width versus thickness shown separately for sheets with X-points (left) and sheets without X-points (right). 
The sheets with X-points appear to have a stronger correlation between width and thickness, with a power law index near 4.0 it 1.0 (similar 
to the length versus thickness fit, Fig. |ll[ |. The width and thickness are not strongly correlated in sheets without X-points, with a possible 
power law index near 2.2 it 1.0. 
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cantly shallower than the fit for the sheets with X-points. 
The difference in cross-correlations of the dimensions be- 
tween sheets with X-points and sheets without X-points 
suggests that reconnection significantly affects the geom- 
etry of the current sheet. 

The scalings of ohmic energy dissipation rate with the 
current sheet length and width are shown in Fig. |13[ 
The energy dissipation rate has a strong correlation with 
both these quantities. The scaling with length is approx- 
imately quadratic, £ oc Z^-^**^-^, and the scaling for the 
width is similar, £ oc ^2.o±o.2_ current sheets, regard- 
less of whether or not they contain an X-point, exhibit 
this scaling. We also find that the energy dissipation 
rate is nearly proportional to the product of length and 
widt h, w ith a power law fit of f oc {IS,)^'^^^'^ , shown in 
Fig. [13] This suggests that neither thickness nor peak 
current density are important parameters in determin- 
ing the total energy dissipation rate from a given current 
sheet. Indeed, although we find a correlation between 
the thickness and energy dissipation rate, as shown in 
Fig. 15 the correlation is relatively weak when compared 
to omer correlations. This could also be expected be- 
cause peak current density and thickness have a smaller 
spread of values than length and width. 

We now consider the scaling of length I with peak cur- 
rent density jmax, separately for current sheets with and 
without X-points. The scatter plot for this is shown in 
Fig. |16[ There is no significant correlation between jmax 
and I for current sheets without X-points. However, there 
is a moderate correlation between the two quantities for 
current sheets with X-points. A power law fit gives a re- 
lationship of I oc Jniax I* interesting that the scaling 
of energy dissipation rate with length is not significantly 
different for current sheets with and without X-points, 
but the scaling of peak current density with length does 
depend on the populations. This may indicate that the 
energy dissipation rate is essentially a geometric quan- 
tity, while the peak current density is more strongly in- 
fluenced by dynamical effects such as reconnection and 
flow into an X-point. 

We find that there is a correlation between the in- 
tensity of the current sheet (measured by peak current 
density ) an d the volume of the current sheet, as shown 
in Fig. 17 In principle, this correlation must be taken 



into account when associating reconnection with current 
sheets that contain X-points, since a current sheet with 
a larger volume has a higher probability of containing 
an X-point purely by chance. This effect in itself could 
produce the X-point asymmetry in the distributions of 
£ or jinax- To address this, we determine whether the 
probability of randomly containing an X-point is impor- 
tant by using the following method. We redo the analy- 
sis with uniformly scattered random points replacing the 
actual locations of X-points. Then we classify the cur- 
rent sheets by whether or not they contain these random 
points. If current sheets with random points reproduce 
the same asymmetry in the energy dissipation distribu- 
tion as observed for current sheets with X-points, then we 
can conclude that the current sheets are not correlated 
with X-points. Upon performing this test, we found that 
the current sheets with random points cannot account 
for the distributions observed by current sheets with X- 
points. Therefore, strong current sheets are indeed most 
likely to form at the locations of X-points, although this 



does not change the fact that a significant number of 
intense current sheets exist without X-points. 

4.2. Comparison to Sweet-Parker model 

The final part of our anal ysis is a test of the appli- 
cability of the Sweet-Parker ( |Sweet||1958"{ |Parker|[l957 l 
model of reconnection, which predicts the reconnection 
rate and current sheet thickness at a site of magnetic re- 
connection, to X-point-containing intense current sheets 
produced in MHD turbulence. Although other resistive- 
MHD reconnection models exist, we consider the Sweet- 
Parker model for the following reasons. First, numerical 
simulations of rec onnection in resistive MHD have shown 
that Pctschek's (jPetschek] [1964) model fails while the 
Sweet-Parker model provides a robust solution, at least 
in regime of modest Lundq uist numbers (Biskamp^^l986 
Uzdensky fc Kulsrud||2000 ). Second, while it is now wel 



established that Sweet-Parker current sheets with high 
Lundquist numbers (S > Sc ^ 10^) become unstable 



to the plasmoid (secondary tearing) instability (Loureir^ 
et al.|[20 07), which c ompletely di srupts th em (Samtaney 
et ai.||2 009 Bhatta charjee et al. 2009; H uang & Bhat- 
tacharjee„2010i ,Loureiro et al. 2012) , the current sheets 



m our simulations are never long enough for their cor- 
responding Lundquist numbers to exceed the instability 
threshold Sc- Therefore, plasmoid-dominated r econnec- 
tion is not relevant in our simulations. Finally, |Servidio] 
et al.| (2009 ) have shown that the Sweet- Parker niodel can 
locally describe reconnecting current sheets in 2D MHD 
turbulence, making it reasonable to hypothesize that 3D 

current sheets behave simil arly. 

The Sweet-Parker model ( |Sweet|1958] |Parker|195"7| as- 
sumes a 2D current sheet of roughly uniform (or slowly 
varying) current density in a high-aspect-ratio, approxi- 
mately rectangular region called the reconnection layer. 
The reconnection rate and the reconnection layer thick- 
ness can then be derived entirely fr om con servation laws 
(for a derivation, see, for example, Biskamp 2003). The 
classic Sweet-Parker model predicts that tor symmet- 
ric reconnection with a corresponding current sheet of 
width ^ and an upstream reconnecting magnetic field of 
strength B, the layer thickness should be 



A 



SP oc 



1/2 



(19) 



where 5* = £,B/ri is the Lundquist number for this layer 
and B is taken in units of the Alfven velocity. A slightly 
modified form of this equation is obtained by estimating 
B in terms of the peak current density in the current 
sheet, jmax « S/A, which gives 



Asp oc 



J71 



1/3 



(20) 



We test the theoretical prediction of the Sweet-Parker 
model for our population of X-point-containing current 
sheets. We use measurements of width, ^, and the mean 
increment in upstream magnetic field, B — \Bi — i?2|/2, 
where the fields Bi and B2 are the in-plane cross-section- 
parallel components of the magnetic field at the two 
edges of the current sheet, as explained in Subsection |3.2| 
(see Eq. 18). We use these measurements to compute 
the expected thickness Asp from Eq. [19] We then com- 
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Figure 13. For all current sheets, there is a strong correlation between the energy dissipation rate and length (left), as well as for the 
energy dissipation rate and width (right). The relationships are approximately power laws, with indices 2.1it0.2 and 2.0it0.2 respectively. 

the center of the current s heet. For upstream paraUel 
magnetic fields Bi and B2, |Cassak fc Shay, ( ,2007| ) find 




Acs « 



1/2 



B2 



1/2 



B2 

Bi 



1/2 



(21) 
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Figure 14. The energy dissipation rate and product of width 
and length are nearly proportional, with a good fit being a power 
law of index 1.1 it 0.1. 

pare Xsp to the actual (directly measured) thickness A 
for the same current sheet. A scatter plot for this com- 
parison, using current sheets fro m th e Case 2 simulation, 
is shown in the left panel of Fig. 18 We observe a large 
amount of scatter around the expected trend. Howe ver, 
if we compute the Sweet-Parker thickness from Eq. [20 
using the measured peak current density jmax insteada 
the magnetic field jump, then we find a better agreement 
and a tighter correlation between the measurements and 
the predictions, shown in the right panel of Fig. [18] This 
suggests that the Sweet-Parker model actually provides a 
reasonable physical picture for what happens in these in- 
tense current sheets. In the other two simulations (Cases 
1 and 3) the agreement is not as strong, which suggests 
that the current sheets in those cases may not be suffi- 
ciently resolved. 

We note, however, that the original Sweet-Parker the- 
ory is intended for describing only the so-called symmet- 
ric reconnection case, Bi — — i?2- Current sheets that are 
formed in MHD turbulence have no a priori reasons to 
be symmetric, however. Therefore, we also consider the 
generalization of the Sweet-Parker mode l to asymmetric 
(Bi ^ —B2) current sheets, developed by Cassak & Shay 
(2007). In this case, the X-point does not coincide with 



We test this prediction in a similar way to the symmetric 
Sweet-Parker model. For the upstream values of Bi and 
B2 , we take the measured fields and subtract off the same 
component of the magnetic field at the current density 
peak. This removes the effect of the large scale field 
and improves agreement. The results are very similar 
to our comparison with the Sweet -Pa rker model using 
B = {Bi- B2)/2, as shown in Fig.[l9[ 

There are a number of reasons why one may expect a 
certain degree of disagreement and a lack of tight corre- 
lation between the measured thickness and the thickness 
predicted by the Sweet-Parker and Cassak-Shay theories. 
In particular, turbulence distorts the current sheets sig- 
nificantly, resulting in violation of some of the assump- 
tions about the current sheet geometry that go into these 
simplified theoretical models. In addition, these mod- 
els are stationary, whereas the current-sheet structures 
developing in our simulations are transient and highly 
non-stationary. 

5. DISCUSSION AND CONCLUSIONS 

This completes our statistical analysis of dissipative 
currents sheets in driven, incompressible, reduced-MHD 
turbulence. We have arrived at a detailed picture of the 
typical current sheets that arise spontaneously from non- 
linear interactions. This analysis has led to the following 
conclusions. 

1) In MHD turbulence, a significant fraction of all 
ohmic dissipation takes place in highly intermittent, rel- 
atively extended structures — in current sheets that are 
larger and more intense than the typical dissipative scale 
naively predicted by MHD turbulence theory. These 
structures are relatively rare and sparse but account for 
a significant fraction of all dissipation. For example, in 
Case 2, more than 25% of all dissipation takes place in 
current sheets with central current peak density higher 
than approximately 8 times the global root-mean-square 
current density; these current sheets occupy less than 1% 
of the volume. 
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Figure 15. For current sheets with X-points, there is a weak correlation between the energy dissipation and thickness (left), fit to a 
power law with index 7.0 it 1.0. For sheets without X-points, the correlation is similarly weak, fit to a power law with indices 4.5 it 1.5. 
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Figure 16. Scatter plots of length versus peak current density, separately for sheets with X-points (left) and without X-points (right). 
Current sheets with X-points have a significant correlation with a power law index 1.2 it 0.2, while there is no strong correlation for current 
sheets without X-points. 
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Figure 17. There is a significant correlation between volume and peak current density for current sheets with X-points (left), although 
there is no strong correlation for sheets without X-points (right). 
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Figure 18. Comparison of measured current sheet tliicknesses to Sweet-Parker theory predicted thicknesses. Iri the left panel, the 
Sweet-Parker thickness is computed using measured width and upstream magnetic field B = (Bi — B2)/2 (see Eg. |19| . In this case, there 
is a large amount of scatter about th e expected values. In the right panel, the Sweet-Parker thickness is computed using measured width 
and peak current density (see Eq.|20[l. In this case, there is better agreement between the measurements and predictions. 
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Figure 19. Comparison of measured current sheet thicknesses 
to the predicted thickness for a symmetric current sheets in the 
Cassak-Shay model (see Eq. |21| l. The results are similar to the 
Sweet-Parker t heor y using an upstream magnetic field of B = (Si — 
B2)/2 (see Fig.jTs}. 

2) Even though there is a natural tendency to automat- 
ically associate current sheets with reconnection, we find 
that this is not necessarily correct. In particular, we find 
that about 55% of the intense current sheets, accounting 
for approximately 16% of the overall energy dissipation 
from current sheets, do not contain magnetic X-points 
(whose presence is seen as a hallmark of reconnection). 
Likewise, we find that there is a large number of X-points 
in the simulation domain at any given time that are not 
associated with strong dissipation sites. Thus, one must 
be careful in postulating a direct correspondence between 
reconnection sites and energy dissipation sites in MHD 
turbulence. 

3) At the same time, however, we find that the most in- 
tense current sheets tend to be reconnecting ones. Thus, 
whereas only about half of the current sheets are found to 
contain X-points, their fraction increases as one considers 
more intense dissipation events. For example, the prob- 
ability of a current sheet to contain an X-point increases 
to 90% when its corresponding peak current density be- 
comes roughly half of the global maximum. Furthermore, 
these current sheets appear to exhibit different statistical 
properties in general; for example, their properties are 
more strongly correlated with thickness and peak cur- 
rent density. These current sheets appear to conform, on 
average, to the Sweet-Parker reconnection theory. 

5.1. Implications 

Our results have potential implications for turbulent 
energy dissipation and plasma heating in many space and 
astrophysical environments, including the solar wind, ac- 
cretion disks and flows, the interstellar medium (ISM), 
and the hot coronal gas in galaxy clusters. 

In particular, these findings may strongly influence our 
understanding of thermodynamics of space and astro- 
physical plasmas where the temperature structure is con- 
trolled by a balance between heating by turbulent energy 
dissipation and radiative cooling, such as accretion disks 
and the ISM. Quite often one treats these environments 



as having a well-deflned quasi-uniform (e.g., slowly vary- 
ing in space) temperature, determined by balancing the 
volume-average dissipation rate with the average radia- 
tive cooling rate. However, this assumption is valid only 
if either the heating is uniform or if the heat released at 
localized dissipation sites is able to distribute itself more 
or less uniformly throughout the system, e.g., via ther- 
mal conduction, on time scales faster than the radiative 
cooling time. If, as we in fact found in the present study, 
turbulent heating is highly localized and intermittent, 
then, depending on the radiative cooling function and 
on the efficiency of thermal conduction, it is in principle 
possible that the plasma does not have a quasi-uniform, 
smooth temperature field. Instead, the system is charac- 
terized in this case by an essentially inhomogeneous and 
intermittent thermal structure, with constantly appear- 
ing and disappearing relatively sparse hot regions that 
are quickly heated by dissipation in current (and vortic- 
ity) sheets and then quickly cooled by radiation cooling, 
before they can share their heat with nearby regions via 
thermal conduction. 

The present work also has potentially important impli- 
cations for situations in which the most intense current 
layers, responsible for a significant fraction of the overall 
dissipation, undergo a transition to fast coUisionless re- 
connection. This may enhance the intermittent character 
of energy dissipation even further. 

Our results may also have implications for the solar 
wind. Magnetic discontinuities, characterized by rapid 
spatial variation in the magnetic field vector, have been 
measured in-situ by spacecraft in the solar wind. The 
observed probability distribution of discontinuities in the 
solar wind have bee n closely reproduced by simulations 
of MH D turbulence ( [Greco et aL||2010[ [Zhdankin et al. 
20121. In our present study, we found that reconnectmg 



current sheets are indeed marked by a strong discontinu- 
ity in the magnetic field, and are reasonably consistent 
with Sweet-Parker current sheets. This may suggest that 
the solar wind contains intermittent current sheets that 
dissipate a significant fraction of the magnetic energy. If 
a more direct association between magnetic discontinu- 
ities (possibly measured by using several different meth- 
ods) and properties of the corresponding current sheet 
can be established, then the solar wind turbulence and 
energy dissipation can be understood in greater detail. 

5.2. Comparison to similar studies 

In this subsection, we compare our results to simi- 
lar statistical studies of intermittent structures in MHD 
turbulence. Our study was preceded b y and partially 
motivated b y the statistical analysis of [Servidio et al. 
(2009 2010), performed on magnetic reconnection sites 
in simulations of 2D MHD turbulence. Their algorithm 
is, broadly, as follows: the first step is to detect X-points, 
and the second step is to study the surrounding dissipa- 
tion region (reverse of our procedure). They estimate 
the two dimensions of the diffusion region by using the 
eigenvectors and eigenvalues of the Hessian of magnetic 
potential. Hence, it is a local estimate that may miss in- 
formation about the current sheet shape away from the 
X-point. In order to be accurate, this approach requires 
a very well resolved current density peak near the X- 
point. Our technique for measuring dimensions is more 
direct and is applicable for a wider variety of current- 
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sheet shap es (for example , highly curv ed sheets). Al- 



though the Servidio et al. 



( |2009l["2010[ ) algorithm is an 
excellent approach tor studying tne reconnection process 
itself, an important limitation is that it ignores current 
sheets not directly associated with an X-point. Our pro- 
cedure extends this idea, making it applicable to dissipa- 
tion sites not involved with reconnection. 



Servidio et al. (2009 20101 find that if the reconnec- 
tion rate at an A-pomt is greater than a fixed threshold, 
then generically the X-point is surrounded by a current 
sheet. For X-points with a low reconnection rate, they 
find no clear scaling properties of the diffusion region pa- 
rameters, suggesting that current sheets are not discern- 
able. For X-points with a high reconnection rat e , they 
find scalings consistent with the Cassak & Shay (20071 
predictions for asymmetric Sweet-Parker reconnection. 
The upstream magnetic fields Bi and B2 and the cur- 
rent sheet length (the width in our 3D terminology) are 
used to test the theory, similar to our analysis. They ob- 
tain a stronger agreement with the Sweet-Parker model 
than we do, which may be due to using 2D simulations. 
Another advantage of their study is that the dissipation 
range in 2D is significantly better resolved than is possi- 
ble in 3D simulations. 



Uritsky et al. (2010) performed a statistical analy- 
sis that more directly studies intermittent structures in 
MHD turbulence. Although the analysis in our present 
paper is in many respects similar to that in their pa- 
per, there are also several important differences. First, 
whereas their paper considers decaying MHD turbulence, 
our present study focuses on driven MHD turbulence. 
Our approach then has the advantage of being able to 
study turbulence in a statistical steady state; our simu- 
lations can in principle run indefinitely, which allows us 
to collect the statistical data using many snapshots and 
thus obtain a much larger statistical sample. The second 
difference is that, whereas our study focuses on reduced- 
MHD turbulence, valid in the presence of a strong guide 
field, their investigation does not have this limitation; it 
uses the full incompressible MHD equations and does not 
assume a strong guide field. In this sense, methodolog- 
ically, their procedures are more universally applicable 
than ours. Thirdly, their cluster algorithm detects not 
only current sheets but also vorticity structures, which 
enables them to study the statistical properties of viscous 
dissipation, in addition to ohmic dissipation of interest to 
us here. In this sense their analysis is more comprehen- 
sive, although one does not expect significant differences 
between the two dissipation mechanisms in the case of 
MHD turbulence, as their study in fact confirms. Adding 
viscous dissipation diagnostics to our study is a straight- 
forward generalization of our algorithm and is something 
we plan to implement in the near future. 

Anot her important differ ence in terms of diagnostics 
is that Uritsky et al. (2010) do not report the statistics 
of the peak current density and do not explicitly men- 
tion the value of the current density threshold used in 
their study; thus, it is not clear to what extent their con- 
clusions apply to the most intense dissipative structures 
(which constitute the main focus of our present study) . 
Also, whereas their algorithm measures several geomet- 
ric properties of the dissipative structures (such as the 
length), the current sheet thicknesses are not measured 
directly. Instead, they are estimated as volume divided 



by surface area, which is, however, a valid approximation 
as long as th e structures are sheet-like. Finally, |Uritsky| 
et al. (2010) do not explore or discuss any association 
of turbulent dissipative structures with magnetic recon- 
nection sites, which is the main goal of our paper. In 
particular, they do not distinguish between reconnecting 
and non-reconnecting current sheets, as marked by the 
presence or absence of magnetic X-points. 

Despite these differences, there are also many i mpor- 
tant simil arities and points of comparison between |Urit-| 
sky et al. (2010) and our study. We both find and ana- 
lyze power-law distributions and scaling relationships for 
many of the measured current-sheet properties. For some 
of them we can make a direct comparison and we find 
that our results are generally in reasonable agreement. 
For example, in our study we find that the current-sheet 
length and width are distributed with a power-law of in- 
dex of about -2.5, which is not too different from their 
value of —tl — —2.2 ± 0.22 for their Run I (see their Ta- 
ble II). Likewise, the ohmic dissipation rate integrated 
over a current sheet has a power-law distribution with 
an index of about -1.8 in our study, which is comparable 
to —1.44 ± 0.06 in their study. We find that the ohmic 
dissipation rate scales with length as a power law with in- 
dex 2.1, while they find an index 2.44±0.10. In addition, 
we both find that the distribution of current-sheet thick- 
nesses A is significantly different from that of most other 
quantities — it declines very steeply at large A, which 
may be more consistent with some sort of exponential 
cut-off rather than with a genuine power law. This re- 
flects the fact that current sheet thicknesses vary across 
a much narrower dynamic range of values and hence do 
not correlate very strongly with, e.g., the current-sheet 
length. 

6. FUTURE DIRECTIONS 

The study reported in this paper should be seen as 
only a first step in a much broader research program 
devoted to understanding the interaction between tur- 
bulence and reconnection in magnetized plasmas. The 
statistical analysis methods we developed can be further 
improved and extended in the following ways: 

(1) Supplementing our existing magnetic-field/current 
diagnostics with an analogous velocity /vorticity diagnos- 
tics, which will enable us to study the statistical proper- 
ties of viscous dissipation in addition to ohmic dissipa- 
tion. 

(2) Adding a temporal dimension to our study, i.e., 
tracking current sheets not only in space but also in 
time, determining their lifetime in addition to their spa- 
tial extent. This should allow us to treat the dissipative 
spatio-temporal structures as events (fiares or flashes) 
rather than just spatial structures and hence investigate 
the statistics of their total (time-integrated) energy re- 
lease, lifetimes, and related quantities. By tracking the 
evolution of individual current sheets in time, we should 
be able to determine whether they coalesce and decay, 
and to assess the validity of a stationary model for re- 
connection. Our present simulation data sets are not re- 
solved well enough in time to perform such an investiga- 
tion. In principle, however, this can be done in the future 
by either having a very high cadence of saving numerical 
snapshots to the disk or by embedding the current-sheet 
diagnostic procedures directly into an MHD code itself 
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and thus doing the analysis simuhaneously with the sim- 
ulation. 

(3) The present statistical analysis can also be im- 
proved upon by using data from better resolved simu- 
lations, which should increase the accuracy of measure- 
ments and reveal stronger trends over a wider sampling 
range. In particular, this could strengthen the corre- 
lations associated with peak current density and thick- 
ness. Indeed, it is known that accurate measurements of 
higher-order statistics require a very well resolved dissi- 
pation range, which is particularly important for numer- 
ical studies of reconnection (see, e.g., |Wan et al.||2010[ ). 
It is difficult to determine whether the current sheets m 
our present study are adequately resolved until improved 
simulations become available. Also, we plan to perform 
a sequence of studies with several different values of the 
resistivity and viscosity. This will enable to see which 
statistical properties investigated in the present paper 
exhibit universal behavior, i.e., are independent of the 
magnetic Reynolds number, and also to see whether and 
how they depend on the magnetic Prandtl number. 

Finally, we expect that the diagnostic methods and 
algorithms developed in the present study will find use- 
ful applications (perhaps with modest modifications) in 
various specific areas of modern heliospheric physics and 
astrophysics, such as: 

- other types of MHD turbulence — most notably, turbu- 
lence driven by the magneto-rotational instability (MRI) 
in astrophysical accretion disks; 

- chaotic magnetic fields in magnetically-dominated en- 
vironments, e.g., in the solar corona in the context of the 
coronal heating problem (and also in coronae of accretion 
disks, e.g., in black-hole systems) ; 

- relativistic turbulence ( Zrake fc MacFadyen||2013 l; 

- non-MHD, purely hydrodynamic turbulence — namely, 
quasi-lD dissipative structures (vortices) in incom- 
pressible turbulence and quasi-2D dissipative structures 
(shocks) in comp ressible superso nic turbulence, e.g., in 



upe 

molecular clouds (|Boldyrev|[2002[) 
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